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INTRODUCTION 


The RBQQ rich zone will be characterized by conditions that exacerbate thermal radiation 
loads on the liners. High soot densities and temperatures in the rich zone may lead to unacceptably 
high radiative fluxes which compromise liner durability. The task of predicting these radiative 
fluxes is complicated not just by the lack of good models for soot formation in jet fuel, but also by 
flow turbulence, which is known to lead to enhancement of radiation. Basing radiative calculations 
on time-averaged CFD temperatures and species densities will generally not be accurate. The 
development of a useful analytical tool for thermal heat transfer prediction thus requires a jet fuel 
kinetics scheme, a soot formation model, and an efficient way of calculating radiative fluxes in 
turbulent environments. The soot model is needed not just for prediction of radiation loads, but 
also for prediction of soot burnout in the quench zone. A further requirement is geometric flexibility 
for the radiation algorithm. 

In response to this problem, a joint UTRC-University of Connecticut theoretical program 
was put in place. The program was based on describing coupled soot formation and radiation in 
turbulent flows using stretched flamelet theory. The University of Connecticut had responsibility 
for Subtask F, entitled Reactive Flow Modelling, and consisting of three parts: development of 
an engineering model of jet fuel kinetics appropriate to diffusive combustion, improvement of the 
standard, linear k-€ turbulence model which is common to many flow and combustion codes, and 
development of a joint pdf methodology for the calculation of mean flow and radiation in a turbulent 
flame. UTRC had responsibility for Subtask G, entitled Flamelet Kinetics and Turbulent Radiation 
Model. This effort was involved with using the model jet fuel kinetics mechanism to predict soot 
growth in flamelets at elevated pressure, to incorporate an efficient model for turbulent thermal 
radiation into a discrete transfer radiation code, and to couple the soot growth, flowfield, and 
radiation algorithms. The soot calculations used a recently developed opposed jet code which 
couples the dynamical equations of size-class dependent particle growth with complex chemistry. 


1 



Several of the tasks represent technical firsts; among these are the prediction of soot from a detailed 
jet fuel kinetics mechanism, the inclusion of pressure effects in the soot particle growth equations, 
and the inclusion of the efficient turbulent radiation algorithm in a combustor code. A schematic 
overview of the main technical tasks and how they are coupled to provide predictions of radiative 
fluxes is shown in 
the two Subtasks, 


the accompanying figure, followed by detailed descriptions of the work done in 
the soot growth/radiation computer program, and sample calculations. 
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Subtask F. Reactive Flow Modeling 


The reactive flow modeling task consisted of three parts: development of a jet fuel kinetics 
model appropriate to diffusive combustion, improvement of the standard, linear k-e turbulence 
model which is common to many flow and combustion codes, and development of a joint pdf 
methodology for the calculation of mean flow and radiation in a turbulent flame. The underlying 
combustion model is based on flamelet theory and is described below. 

FLAMELET MODELING 

Turbulent combustion modeling necessitates dealing with the description of reaction rates, 
directly or indirectly. Advanced combustion models seek to circumvent the problem of solving 
balance equations and associated closure problems involving scalar correlations which appear in a 
conventional Reynolds decomposition of the species conservation equations. 

Of the more advanced methods flamelet and PDF methods are the most promising. PDF 
methods potentially have greater generality, but are mathematically more complex than flamelet 
models and presently are limited to very simple chemistry. The latter limitation precludes the 
modeling of soot formation from large hydrocarbon molecules. Modeling of turbulence/radiation 
effects using PDF methods also is expected to be especially cumbersome. On the other hand, 
flamelet modeling, the approach followed in this work, easily accommodates complex chemistry and 
radiation effects. Although in certain applications the principal constraint in flamelet modeling 
may be a physical one, that the scale of the reaction zone need be small relative to the scale of 
turbulence, this is not expected to be a limitation in modeling high pressure, high temperature 
flames as found in gas turbine combustors. 

In flamelet modeling, the combustion zone is treated as an ensemble of folded, laminar-like 
structures which are convected by the turbulence. The thickness of the laminar-like reaction zone is 
a function of the strain rate as reflected in the molecular (scalar) dissipation, a quantity analogous 
to the viscous dissipation. Both the mean scalar and viscous dissipation are computed in the mean 
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flow calculation. 


A highly significant advantage of the flamelet approach is that kinetic calculations are not 
carried out in the main flow calculation. Thus unencumbered, the problem is reduced to computa- 
tion of mixing with variable density. Luminous and non-luminous radiation, intensity combustion 
products, and soot concentrations are derived from post-processing of the main flow data using 
analytical expressions for these properties derived from a laminar flamelet calculation. 

JET FUEL KINETICS 

The objective of this task was to develop a kinetics scheme which provides realistic estimates 
of flamelet temperature and flame products concentrations necessary for the prediction of both 
luminous and non-luminous radiation. The prediction of intermediates, particularly acetylene and 
aromatics, was considered essential for the description of soot formation. To keep this task man- 
ageable under the allotted time it was necessary to do two things. The first was to decide upon 
the composition of the simplest model fuel which might mimic the behavior of jet fuel, and the 
second was to devise a kinetics scheme for the model fuel. Both pyrolysis and oxidation kinetics 
were considered. Regarding the latter, simple, global kinetic schemes, as described in the literature, 
were not considered adequate for two reasons, the lack of description of intermediates necessary 
to describe the formation of soot, and the limited range of pressure, temperature and composition 
over which global schemes are valid. Furthermore, diffusive effects, not reflected in global schemes, 
are important in practical devices. The development of “jet fuel kinetics” therefore necessitated 
considering both complex pyrolysis and oxidation kinetics keeping in mind that combustion would 
occur primarily in a diffusive (as opposed to premixed) mode. This latter assumption allowed a 
major simplification of the kinetics scheme. The kinetics scheme and its inclusion in the flamelet 
calculation is discussed below. 
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MODEL FUEL SELECTION 


One important consideration in devising a model fuel is that the fuel should mimic the com- 
bustion of jet fuel in the formation of soot. The appropriate physical model of soot formation is 
that soot derives from the inception and subsequent growth of soot nuclei originating in fuel-rich 
zones. The inception process consists of the formation of polycyclic aromatics (PAH) under fuel-rich 
conditions followed by the growth of PAH from the continuous addition of acetylene to the PAH. 
Eventually a solid phase is formed which also grows grow by acetylene addition. Soot kinetics are 
discussed in Subtask G, where a detailed description of the soot formation model used in this work 
is given. 

Since jet fuel contains about 20% aromatics, incepting species are abundant initially, and the 
inclusion of an aromatics component in the model fuel is essential. Acetylene is a product of the 
pyrolysis of the alkane constituents, the most abundant family of compounds in jet fuel. Although 
the alkane fraction consists of hundreds of compounds, individual alkanes will pyrolyze in a similar 
manner to yield ethylene, methane, hydrogen, and most significantly, acetylene. Thus, a simple 
model fuel would consist of two classes of compounds, alkanes and aromatics, and the simplest 
model fuel would contain a single representative alkane and a single representative aromatic. 

Detailed chemical analyses provided by Southern Petroleum Laboratories and Pratt and Whit- 
ney Aircraft were the principal sources considered in devising the model fuel. Total saturates and 

aromatics and the most abundant carbon number for each class are shown in Table 1. 

Table 1 



P&W 

SPL 

Average Carbon Number 

Most Abundant 

Saturates 

79.1 

76.6 

9.98 

Cio 

Aromatics 

20.5 

19.1 

9.12 

Cg, Cio 


N-decane and trimethylbenzene were selected as the model fuel components on the basis of 
average carbon number and the detailed analysis (85% n-decane, 15% trimethylbenzene). Contri- 
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butions to soot from other relatively abundant compounds, such as indans, tetralin, and napthalenes 
were not considered as these could be considered as additive incepting species (see soot model). Cy- 
cloalkanes were not considered separately. 

The overall kinetic scheme is summarized in Fig. 1. The reaction equation set is presented in 
Appendix A. This scheme was incorporated into the laminar flamelet calculation which provided 
flamelet temperature and gaseous species distributions. In a second step these data were used as 
the starting point for the calculation of soot via the UTRC soot model. 

For n-decane, it was assumed that appreciable heating of the fuel in a diffusion flame occurs 
prior to exposure to oxygen. Thus a reasonable simplification would allow that substantial decane 
pyrolysis occurs before significant oxidation. The process was viewed as one in which decane served 
as a source of pyrolysis products, most of which were subsequently oxidized and a small fraction 
of which (C 2 H 2 ) participated in the formation of soot. Oxidation was confined to Cl and C2 
species. The n-decane kinetics are modeled through a step-wise pyrolysis scheme beginning with 
the formation of decyl radical and proceeding to the formation of smaller saturated and unsatu- 
rated molecules whose oxidation is modeled comprehensively. An abbreviated benzene formation/ 
pyrolysis/oxidation scheme was used. The approach is opposite to that taken in descibing decane 
kinetics. It was assumed that aromatic nuclei were thermally stable. Aromatics kinetics was mod- 
eled using a global oxidation scheme for 1,2,4 trimethyl benzene (TMB), the most abundant jet fuel 
aromatic constituent (Ref. 1). The mechanism describes the stepwise oxidation of mono-aromatic 
intermediates and decomposition to benzene (see Appendix A). 

LAMINAR COUNTERFLOW FLAME (FLAMELET) CALCULATIONS 

The kinetics scheme is incorporated into a laminar, counterflow, diffusion flame calculation 
which serves as the basis of the turbulent combustion model. The flamelet calculation is done 
for a counterflow diffusion flame with appropriate temperature, pressure, and mass flux boundary 
conditions for the fuel and oxidant streams. Mixture fraction is computed from the species profiles, 


7 




products 


Figure 1. Model jet fuel kinetics schematic. 
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and the relevant species, temperature and density functions are expressed as polynomials in mixture 
fraction. Individual flamelet calculations are parameterized by the scalar dissipation rate which is 
a function of the reactant mass fluxes. Successful laminar flame solutions have been obtained 
with the model fuel. Substantial concentrations of benzene and acetylene are indicated. Typical 
flamelet profiles are shown in Fig. 2. Inspection of the data computed at typical engine operating 
conditions, indicate that the the heat release zone is very thin, about 0.3mm, thus supporting the 
flamelet model. 

NON-LINEAR k-e TURBULENCE MODEL 

A non-linear k-e turbulence model, based upon the work of Speziale (Refs. 2, 3), was incorpo- 
rated into the TEACH code. The non-linear model provides a significant improvement in predictive 
capability over the standard k-e model without significant additional computational burden. No 
additional equations are introduced into the analysis as with other advanced models such as dif- 
ferential Reynolds stress (DSM) algebraic stress (ASM) models. Notable successes of the model 
include improved prediction of normal stresses in channel flows, prediction of secondary flows in 
non-circular ducts (inherently impossible with the standard linear model), and improved prediction 
of recirculation zone length behind a rearward facing step. 

The standard, linear k-e model, assumes a linear relationship between stress and mean vorticity 
rendering it inherently unable to describe secondary and other flows with anisotropic normal stresses. 
For curved or swirling flows, correction terms are invoked in the dissipation equation and eddy 
viscosity formulation. Often these solutions are found to be problem dependent ,and thus, limited 
in generality. More advanced, (and also more complicated) approaches such as ASM and DSM 
may offer greater potential for complex flows, but this has not been established, and the additional 
computational burden often may not be justified. The non-linear stress model extends the validity 
of linear stress models by allowing for normal stress anisotropy. Two additional quadratic terms are 
added to the momentum stress expression. The additional terms are subject to several mathematical 
constraints, the principal restraint being that of frame indifference (Refs. 2, 3). This requires that 
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Figure 2. Species profiles in model jet fuel flamelet. 
Fuel: 15% trimethylbenzene, 85% decane 


the tensor form of the non-linear terms not change with change of reference frame, i.e. the form is 
the same for both inertial and non-inertial reference frames. The expression is written below. 

2 k 3 

nj = -r/Mj + 2Cfip — Dij 

O t 

^3 r 

+CDp(2C/i) 2 --^ j^Dim D m j — l/3D mn Dmn^ij (l) 

+v • VDij - -^^-Dkj - ~~ Dki 
dxpt ^Xk 

+ 'g'[ v ’ - 2 ^k Dmk l 

where Dij is the mean deformation tensor. 

The constant C D is assumed to have the value 1.68. The first two terms on the RHS represent 
the usual, linear formulation for the momentum stress; the remaining terms comprise the non- 
linear contribution. The expression above was rewritten for generalized coordinates and applied to 
a two-dimensional cylindrical system for application to the TEACH code. The principal difficulty 
in implementation is that a large number of quadratic terms are introduced when non-Cartesian 
coordinates are used. 

Predictions of the linear and non-linear models are compared for a model fuel combustion in 
an axisymmetric, constant radius, dump combustor configuration (Fig. 3). Profiles of temperature, 
velocity, and soot volume fraction (soot volume fraction shown as a function of mixture fraction) , 
computed by the joint pdf method, are compared in Figs. 4-7 for the following conditions: pressure, 
10 atm, air temperature 917K, fuel temperature 478K, and equivalence ratio, 0.5. The predictions 
of the linear and non-linear models are seen to differ considerably, indicating that the net radiative 
loss and other flame properties may be sensitive to the choice of turbulence model. Also shown is 
the soot volume fraction derived from laminar flamelet data for a single value of strain rate. Since 
the absolute level of soot volume fraction is strain rate dependent the most significant finding in this 
comparison is that one effect of turbulence is to shift the peak volume fraction to leaner mixtures. 
More numerical examples will be given in the Subtask G discussion. 
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Figure 3. Model dump combustor. 
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Figure 4. Comparison of centerline temperature profiles in TEACH code simulation. 
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R=0.0 m 



Figure 5. Comparison of centerline velocity profiles in TEACH code simulation. 
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Figure 6. Joint pdf and laminar flamelet state relationships compared. 
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Figure 7. Joint pdf and laminar flamelet state relationships compared. 


16 


MODELING THE JOINT PDF 


Mean flame properties are derived from underlying flamelet properties by statistical methods. 
The random variables involved in the averaging process are mixture fraction, z, scalar dissipation 
rate, x> and a characteristic value of scalar dissipation rate which identifies a particular flamelet. The 
value at stoichiometric, x»t> is usually chosen, x and z are assumed to be statistically independent 
(Ref. 4). Experimental evidence indicates that z and x are well represented by beta and lognormal 
probability density functions, P(z;ai), and P(x; A), respectively. a\ and j3\ are each two parameter 
sets which complete the description of the probability densities. Three of the four parameters 
are described by transport equations which form part of the overall mean flow equation set. The 
remaining parameter (see Appendix B), related to the variance of x> is experimentally derived and is 
assumed constant throughout the flow field. An expresion for the joint probability density function 
is presented below. 


Flamelet properties such as temperature, density, and composition are computed in the counter- 
flow flame calculation subject to imposed reactant flow rate boundary conditions which are related 
to x«t • Individual flamelets are identified by their characteristic scalar dissipation rate Xat- Any 
flamelet property Q(z,x*t) will then have a mean local value given by the following expression, 


roo r 1 

(Q) = / / Q(x.t,z) P(x.t,z) dz dx.t 

Jo Jo 


( 2 ) 


where P(x«t> z) is the joint density of z and x»t- 

The joint density is given by the expression (Ref. 4), 


P(x,t> z) = f • P x(x»tf) • P(s) 


( 3 ) 


The shape factor f is defined as x/x«t and is found to be nearly independent of the flow rate 
boundary conditions. Thus, f is a function of z only. 
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The marginal density P(x.t), not used explicitly, but used as a check of the distributon function 
for x.t, is obtained by integrating the joint density in z space. Thus, 


P(x..) = f 

Jo 


P(x»t,z) dz 


( 4 ) 


Numerical Integration 

Since x»t is unbounded it was necessary to assign an upper limit to x.t, x.t max, which yielded 
a value of the distribution function close to unity. The distribution function F(x.t) is defined by, 

rx. t 

F{X ' t) = Jo ( 5 ) 

In principal, the distribution function could be evaluated at each point as a function of x.t max 
M< i X«t max chosen to satisfy some pre-determined value of the distribution function. This procedure 
would have been overly cumbersome, however, so an approximate procedure was used. An estimate 
of the standard deviation of x.t was derived from analytical expressions for the moments of x and 
assumptions of mean values of the shape factor f(z). The upper limit of x.t was then assumed to be 
three or four standard deviations above the mean value of x.t- Occasional checks of the distribution 
function showed that an upper limit so defined was sufficiently large. Further details are given in 
Appendix B. 
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Subtask G. Soot Kinetics and Turbulent Radiation Model 


This task, closely coupled to Subtask F, was comprised of the following: 

1. Modification of a discrete transfer radiation code to provide for boundary-fitted coordinates 
in axisymmetric annular combustor geometries, and the inclusion of realistic, wavelength-dependent 
gas and soot radiative coefficients. 

2. The inclusion of an efficient algorithm for turbulence effects on gas and soot radiation in 
the above code. 

3. The prediction of soot growth in the model jet fuel fiamelets at elevated pressure. 

4. Coupling the output of the CFD solver and the sooting flamelet calculations to the radiation 
code. 

TURBULENT RADIATION 

Turbulent fluctuations enhance time-averaged radiation from flames relative to predictions 
based on time averaged flame properties, and the effect can be very large (Refs. 5-8). An analytic 
treatment of turbulence effects on monochromatic radiation was provided by Kabashnikov and Kmit 
(Ref. 5) for the Wien spectral regime and an assumed linear variation of absorption coefficient with 
gas temperature. Subsequent analysis of the effect in combustion has been mainly numerical in 
nature. The “Monte-Carlo” modelling approach of (Refs. 6,7) divides optical paths into a number 
of homogeneous, statistically independent elements with dimensions corresponding to the turbulence 
integral scale, and sets up possible instantaneous realizations of optical path. This has been done by 
randomly sampling the fuel mixture fraction distribution function within each homogeneous element, 
with the underlying pdf parameters derived from a turbulent flow model solution. Assumed state 
relationships between sampled mixture fraction and temperature/radiating species concentrations, 
usually taken from laminar flamelet solutions, complete the scheme. The inhomogeneous path 
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parameters are supplied as input to standard radiation band models, and the radiative intensity 
pdf is built up by performing many trials. This approach of setting up many realizations may 
not be practical in modelling complex combustor geometries with large numbers of grid points, 
however. The complexity of this approach has seemed so daunting that many modellers have been 
forced to neglect the effect in the hope that it is small in cases of interest to them. A simpler 
and faster semi-analytic approach for gas and soot radiation is described here, and found to give 
good agreement with the more cumbersome Monte-Carlo approach. It rests on the decorrelation 
of point- and path-averaged properties. A simpler calculation results in which attenuation-related 
terms are based on time averaged properties, and the local radiant power density is ensemble 
averaged over the fluctuation pdf using efficient numerical quadratures. Only one path integration, 
yielding the time-averaged intensity, is needed for the spectrally-integrated soot emission, and for 
each molecular band. The result is essentially equivalent to Monte-Carlo with a great reduction in 
computation time. The need to perform pdf averaging of the local radiant power density at each 
node point represents little more effort than is ordinarily expended in turbulent flow calculations 
where ensemble-averaged properties are desired. Numerical examples showing the application of 
this theory to a CH4-H2 turbulent diffusion flame and to a research combustor will be presented. 

In the absence of scattering, and neglecting wall effects, line-of-sight monochromatic radiation 
can be represented by the integral 

I = f k(w,s') I b (a/, s') e’£ k( "’*" ) d ’" (6) 

Jo 

where k and I b represent the local absorption coefficient and Planck function, respectively. For a 
fluctuating medium, the ensemble- or time-averaged intensity is represented by 

(I) = f <k( W ,s') I b (a;, s') e~ k( "’*" ) d *"> ds' (7) 

Jo 

Kabashnikov and Kmit suggested that under certain circumstances it suffices to replace the absorp- 
tion coefficient in the attenuation term with its time average value, and to employ a time-averaged 
local power density at each point, e.g. 


20 



( 8 ) 


(I) ~ f <k(w,9') I b (w,s')> e 

Jo 


<k(*,s")> <«»" dg / 


This is equivalent to saying that the time-averaged intensity is due mostly to fluctuations in 
the local emission power density, with the path-dependent exponential attenuation terms averaging 
out to some extent and making less of a contribution. The simplification afforded by this result is 
obvious, provided that something is known about the local turbulent fluctuation probability density. 
The only alternative is the inherently inefficient process of using random number generators to set 
up realizations of optical path, and building up the radiation statistics by performing a radiation 
calculation for each such realization. Hall and Vranos (Ref. 8,9) arrived independently at a similar 
conclusion for spectrally-integrated, wideband gas radiation. They verified their result for turbulent 
diffusion flame radiation by comparisons with such “Monte Carlo” calculations. (In this paper, 
Monte Carlo will be used in a somewhat different sense than it usually is in radiation calculations). 
Kabashnikov and Kmit, and Krebs, et al (Ref. 10) have shown that a condition for the validity of 
this approach is that individual eddies not be optically thick. For soot radiation, this condition will 
usually be satisfied. For the CO 2 4.3p band, there may be violation at high pressures, but usually 
the soot radiation will be dominating. The application to gas and soot radiation in combustors will 
now be discussed. 


GAS RADIATION 

The analysis of turbulent gas radiation has been given by Hall and Vranos using the exponential 
wideband model (Refs. 8, 9, 11). It will be illustrated by application to a turbulent diffusion flame. 
With the usual simplification that the absorption features vary much more rapidly than the Planck 
function, the intensity generated along a line of sight can be represented as 


I 



dAi 

ds 7 


Ib K 


(o) 


s') ds* 


( 9 ) 
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where I b is the Planck function evaluated at band center frequency u>f 0) , and the subscript i denotes 
the i-th active molecular band. Here Aj is the integrated band absorptance, which for combustion 
problems can be well approximated using the high pressure form 


= (Mf/Aw) + Ei (C/Aw) + 7 b) (10) 

C = a p |s — s'| 

where Aw is molecular resonance bandwidth, o is integrated band intensity, and p is the infrared 
active species density. Band overlaps are neglected, and it is understood that a summation over 
all active bands of HjO, CO 2 , and CO will be performed to calculate intensity. Neglecting band 
overlap effects, the subscript can be suppressed with the understanding that a sum over all bands 
will be performed at the end. 

There is much current discussion in the radiative transfer co mmuni ty about the relative merits 
of narrowband, wideband, and line-by-line calculations (Ref. 12). The latter are too time consum- 
ing for practical applications at the present time, and narrowband models are thought to be more 
accurate than the wideband. However, for engineering purposes the computationally efficient wide- 
band models are felt to be acceptable (Ref. 13), particularly since soot radiation will be dominant 
in most cases of interest. As will be seen in the next section, the calculation of soot radiation is 
more nearly exact, requiring none of the approximations that are inherent in the band models. 

For nonhomogeneous optical paths, the Curtis-Godson scaling approximation as given by 
Morizumi and Edwards (Ref. 14) is employed. In this approximation, the band absorptance is 
expressed in terms of scaled parameters as 


A 

Aw 


tn (£/Aw) + Ei (f/Aw) + 7 e 


( 11 ) 


where 



( 12 ) 


£ = Jap ds" 

AtJ = \ f oj p a ds" (13) 

£ -'s' 

Thus, we have, approximately, 

I ~ J* (^) (1 - e'^) pah ds' (14) 

The quantity of interest is the ensemble average of Eq. 14. In the integrand, the factor pa lb is 
point specific, but multiplies an A-derivative factor that involves only path averaged properties, as 
per Eqs. 11-13. Inasmuch as these paths traverse eddies or volume elements that are presumed to 
be statistically independent, one can make the approximation that the two factors in the integrand 
are statistically independent, i.e. 




{pa h) ds' 


(15) 


where ( ) denotes ensemble- or time- aver aging. 

It will be assumed that state relationships giving temperature T and species densities in terms 
of the fuel mixture fraction pertain. These are obtained in laminar flamelet theory from opposed jet 
solutions at a representative value of strain rate; typically the radiating gas species concentrations 
and temperature are not strongly sensitive to strain rate if the strain rate is in the appropriate 
range. Thus, if the probability density p (z, rji) for mixture fraction z is known, where r]\ are the 
known parameters of the pdf, the ensemble average 
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( 16 ) 


(p a I b ) = f p(z; p(z) a(T(z)) Ib(T(m)) dz 
Jo 

can be regarded as a known quantity at each point along the optical path. The crux of the analysis 
lies in showing that the ensemble average of the other function of path- averaged properties is 
adequately represented by evaluation of the function with time-averaged properties, i.e. 


(jr (1 - ^ (l - (17) 

where the terms on the r.h.s. are evaluated on the basis of time- averaged properties. Details are 
given by Hall and Vranos. 

Example calculations are now shown for an atmospheric pressure, CH«-Hj turbulent diffusion 
flame on which extensive diagnostic measurements have been reported in Ref. 15. The flowfield 
was simulated with a standard k-c turbulence model and a parabolic flow solver, providing at each 
spatial node point the mean fuel mixture fraction and its variance. We assume for these example 
calculations that the mixture fraction pdf p(z) is described by the beta density (Ref. 16), 


P(z) 


T(a + b) , 

r(a) r(b) 


a= 7(z) 




b = 7 (l-<z)) (18) 

a+ b = 7 

.. _ (z) (1 - (z» 

5*5 — 1 

where T is the gamma function, and (z) and (z /2 ) are the mixture fraction mean and variance, 
respectively. The state relationship between mixture fraction and temperature, density, and species 
concentrations was assumed to be given by an opposed jet or counterflow flame solution, employing 
a widely used program (Ref. 17). Mixture fraction is here defined as the average of the C- and 
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H-atom mixture fractions, each of which has been normalized to its fuel side value. These solutions 
are characterized by a representative value of the strain rate, roughly the velocity gradient normal 
to the flame structure, and a solution corresponding to a median value is used. When soot is 
included in the flamelet calculations, the joint pdf of z and scalar dissipation (strain rate) must 
be employed, as has been discussed in Subtask F. The k-e based parabolic code gives a solution 
in distances normalized by the inner fuel tube radius (Figure 8). The optical paths were divided 
into segments of length corresponding to the local integral scale. The “Monte Carlo” calculations 
were then performed in a way similar to Refs. 6, 7. Within each independent volume element, a 
random number generator was used to randomly sample the mixture fraction distribution function 
(Ref. 6). From the resulting value of mixture fraction the instantaneous temperature and species 
concentrations were then interpolated from the opposed jet state relationships. The intensity for 
the realization was then calculated from Eqs. 15-17, summing the active bands. It is also possible to 
make a calculation of radiative flux based on time-averaged temperature and density. The quantity 
of interest is the ratio of time- averaged intensity to intensity based on time-averaged properties, or 
the intensity enhancement. 

Table 2 compares the “Monte Carlo” and “analytic” predictions for the time- averaged, line-of- 
sight intensity at two heights above the burner surface. The two heights encompass a significant 
range of optical thickness, as shown by the band center optical depth of the H 2 O 6.3 fx transition. 
As seen, the “analytic” predictions, which are much more efficiently obtained, satisfactorily agree 
with the Monte Carlo predictions. There are minor differences, not shown, in the absolute value 
of predicted mean properties having to do with the much different algorithms for the two types of 
calculations. This gas band discussion has employed wideband models. If the use of narrowband 
models is preferred, there seems little reason to believe that the underlying result would not be 
applicable to them, as well. 


25 



Divide path into 
homogeneous 
segments 



2Ro 


Figure 8. Representation of turbulent diffusion flame showing division of optical path 
for radiation calculations into statistically independent segments. 



Table 2 


Intensity Enhancement Factors - r( (T57 


H (cm) 

“Monte Carlo” 

“Analytic” 

(C/Awje.s^ 

25 

1.317 

1.327 

0.73 

200 

1.194 

1.197 

2.54 


SOOT RADIATION AND DISCRETE TRANSFER ANALYSIS 

The discussion of turbulent soot radiation is illustrated with an application to a realistic com- 
bustor geometry. A flux model for axisymmetric, annular geometries has been developed using the 
discrete transfer model of Lockwood and Shah (Refs. 18, 19). The discrete transfer algorithm has 
been selected because it gives considerable geometric flexibility, as is well known. A boundary-fitted 
coordinate system that uses transfinite interpolation is employed, so that curved inner and outer 
radial boundaries can be handled. Application of the program to a simulated annular combustor 
geometry is shown in Figure 9. The program works its way around the boundaries of the com- 
bustor, at each point P firing rays in all directions into the combustor, and locating the points of 
intersections with the walls. These wall intersection points serve as starting points for line-of-sight 
radiation calculations back to the point P, at which the net radiative flux is calculated. The x 
symbols denote the points at which the rays pass through axial and radial boundaries; for each 
ray, a list of the cells passed through and the length of the ray within each cell is made for the 
radiation calculation. Because medium properties will vary from cell to cell, this allows medium 
inhomogeneities to be accounted for. The grid should be fine enough to resolve the gradients in 
average temperature; if the turbulent radiation algorithm is to be applied, it is important that the 
properties in adjacent cells be statistically independent. The rays which seem to be highly curved 
correspond to paths emanating from the outer radial wall, miss the inner radial wall, and end in the 
outer wall. The radiative flux divergences (from which gas cooling rates can be derived) can also 
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Figure 9. Discrete transfer rays in model, axisymmetric annular combustor. 



be calculated for each cell; to do so requires calculation of the internal radiative intensities. Wall 
fluxes tend to become relatively insensitive to the number of rays at around 32 rays per point per 
quadrant; depending on the geometry and mesh, however, larger numbers of rays may be needed to 
sample certain remote cells for accurate radiative dissipation calculations in these cells. The line-of- 
sight radiation calculations can be performed either with a narrowband radiation model (Ref. 20), 
or with a combination of the wideband gas model discussed in the last section and a quasi-analytic 
soot model which will now be discussed. 

The calculation of spectrally-integrated soot radiation along a homogeneous path can be repre- 
sented in closed form if gas radiation effects are small. Given typical thermal radiation wavelengths, 
soot particles are usually in the small size parameter or Rayleigh range where scattering is negligi- 
ble. To first order, if intracluster multiple scattering effects are small, the soot absorption coefficient 
can be taken to have the form appropriate to Rayleigh spheroids even for clusters, e.g. 

k 8 (u;) = c 9 u;f v (19) 

where f v is the particle volume fraction, w in units cm” 1 , and the constant c s is related to the 
complex soot index of refraction. It will be convenient to ignore both the frequency dispersion of 
c s and its temperature dependence. If the Planck function is represented as 

oo 

Ife(T) = ciw 3 e - nCau/T (20) 

n=l 

line-of-sight soot radiation can be represented by the double integral over path and frequency 
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the frequency integral can be done first, using 
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giving 
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which is equivalent to, for a homogeneous path, 
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If the optical path is inhomogeneous (Figure 10), the intensity line integral, Equation 21, can 
be represented as a sum of N terms over each of which the path integral can be done analytically, 
as in Equation 24. This gives the inhomogeneous path expression 
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where 

Fi = E f vU)A0) 

j=l 

A wall at Q with temperature T w will add to Equation (25) the expression 


Iw — 3!e w ci ^ 
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If the wall is non- black, and we can assign an effective radiation temperature T r to the flux incident 
on Q, Equation (26) can be replaced by 
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Figure 10. Representation of discrete transfer optical path for 
analytic soot radiation analysis. 



1 


(27) 


Iw = Iw (Eq. 26) + 3!(1 - e w )ci 


(ncj/T r + c,F n ) 4 
as the contribution of the opposite wall to the radiation at P. 


Radiation fluxes calculated using Eqs. 25-27 agree precisely with those using the narrowband 
RADCAL representation (Ref. 20) provided that a consistent value of c» is used (assumed equal 
to 7.0 for c.g.s wavelengths). In the Eqs. 25-27 calculation the frequency summation is done 
analytically as shown; in the narrowband calculation it is done numerically, thus giving rise to 
considerable efficiency enhancement in the former case. A sample comparison is shown in Figure 
11 for the outer annular wall of the model combustor shown in Figure 9. A uniform temperature 
of 2000 K and uniform volume fraction of 10 -6 and cold, black walls have been assumed. The 
agreement is also exact for inhomogeneous medium test cases. 

To calculate soot radiation from turbulent flames, the stochastic analysis requires that fv^) in 
the integrand of Equation 21 be decorrelated from the path integral of volume fraction appearing 
as the exponential argument in the integrand. Doing so (but keeping the correlation of f v (l) with 
the path integral for A(l)) gives 
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with 

<F.) = E (f.(j)>A(i) 

j=l 

where the Tj and the f v (i) are fluctuating quantities. 
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Figure 11. Net radiative flux along outer radial surface in Figure 9 model combustor. 



This has the following limits 


Optically thick 

OO 

<I)^3!cic 2 - 4 (Tj) £ ~ = ~{Tj) (29) 

n= 1 

Optically thin 

oo N 

<I)^4! C1 C 3 - B 22 p E <«r(i)Xf>A(i) (30) 
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In the optically thick limit, the time-averaged flux loses its sensitivity to fluctuations in volume 
fraction, and becomes proportional to the time-averaged value of T 4 at the wall. Wall contributions, 
Eqns. 26, 27 will be evaluated on the basis of average soot volume fraction and the time-averaged 
T r . 


Discrete transfer calculations for an annular combustor shape representative of the RQL rich 
zone are shown in Figures 12 and 13. For the homogeneous temperature and volume fraction values 
chosen, the net radiative fluxes on the inner and outer radial surfaces (assumed to be cold and 
black) are shown. The outer wall sees more radiating gas and has higher predicted fluxes. On 
both surfaces, the contraction is predicted to be an area of high flux because the surface normal is 
exposed to a longer optical path. 

COMBINED GAS AND SOOT RADIATION 

The combined effect of gas and soot radiation is not simply additive, mainly because of ab- 
sorption of gas band radiation by the soot continuum. The gas bands also absorb soot radiation 
to some degree. Each gas band’s radiation is attenuated by the soot absorption coefficient at the 
band center frequency. Thus, Eq. 9 becomes 




where the i summation extends over all molecular bands. 


( 31 ) 
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Figure 12. Discrete transfer analysis of RQL rich zone. 
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Figure 13. Net radiative fluxes on RQL inner and outer radial walls. 


The absorption of soot continuum radiation by the gas bands employs an approximation given 
in Ref. 21. The total soot radiation is Equation 25 plus the expression 


c, £ «< 0) jf ’ d.' f. I b , Awi (32) 

where rj‘* = c,W; 0) / g * f v ds"; ft and Awj are as defined in Equation 7 and 8, and Ibi = Ib(wj (0> ). 

Thus, the total radiative flux when both gas and soot are present is given by Equations 25, 31 
and 32. In turbulent media, the soot and radiation follows the previously discusses treatments, but 
the soot attenuation term in Equation 31 and the entire Equation 32 are based on time-averaged 
properties. This treatment of gas-soot overlap effects is not exact, but its predictions are generally 
close enough to those based on more exact narrowband calculations to be acceptable for engineering 
purposes. 

Soot scattering has been ignored in this analysis beacause typical soot sizes of 0.05 to 0.1 
micron and the wavelengths of interest in thermal radiation analysis produce relatively small size 
parameters (xD/A) that are associated with weak Rayleigh scattering. The large imaginary part 
of the soot index of refraction also dictates that the soot extinction is largely due to absorption. 
Scattering could, however, be included if unusually large soot particles were predicted or observed. 

SOOT DISTRIBUTIONS IN MODEL JET FUEL FLAMELETS 

As discussed, the microflow for the stretched flamelet/pdf approach to turbulence modelling 
is the counterflow flame, depicted in Figure 14. The soot concentrations for the model jet fuel 
flamelets were calculated using an opposed jet solver with complex chemistry that has had added 
to it the dynamical and transport equations of soot spheroid growth (Ref. 22). The modular 
spheroid dynamical model includes inception, coalescence, surface growth, and oxidation, and uses 
the sectional or size bin representation. The model also accounts for soot particle scrubbing of 
growth species and oxidants, and for gas and soot radiation. Before this model could be applied 
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Figure 14. Configuration of opposed jet, laminar diffusion flame. 
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to the present problem, however, a number of code enhancements had to be carried out. The first 
was the creation of a potential flow version more appropriate for flamelet calculations than the 
plug flow assumption used in the original program. Second, since high pressure calculations were 
desired, pressure effects on the soot growth kinetics had to be included. As pressure increases, the 
largest particles move from the free-molecule growth regime (mean free path larger than particle 
size) to the so-called continuum regime. Correction for this effect is non-trivial, and certain of 
the needed code modifications were carred out on this contract. The main correction made was 
to correct the particle coalescence frequencies for continuum effects. This has a pronounced effect 
for pressure levels on the order of ten atmospheres, and it is felt that the corrections made are 
adequate for moderate pressures. To extend the calculations to much higher pressures, however, 
further corrections to the thermophoretic velocities and the surface growth expression would be 
required. 

The model jet fuel kinetics mechanism had 42 species and 123 chemical reactions (Appendix A). 
The soot volume fraction is generally insensitive to the number of size bins, but accurate average soot 
size and number require more bins. A compromise of eight bins was employed in these calculations. 
Volume fraction is the soot size/density parameter of most interest for radiation calculations. The 
major uncertainties in soot loading calculations are associated with the modelling of the nucleation 
and surface growth processes. The field of soot growth modelling is being actively investigated, and 
answers to these questions are being refined. For these calculations, a relatively simple, but widely 
used and highly successful model due to Lindstedt and co-workers (Ref. 23) has been used. While 
developed for methane, it has been applied with good results in many simulations, one of the most 
noteworthy being in co-flow, acetylene-air diffusion flames (Ref. 24). In this model, both inception 
and surface growth are proportional to acetylene concentration, with Arrhenius-type temperature 
dependences. The dependence of inception on acetylene concentration is the most controversial 
aspect of this model. However, models of inception based on small mass PAHs have not yet been 
entirely validated and accepted. The Lindstedt model is felt to be the best available at this time. 
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Another favorable feature of the stretched flamelet approach; wherein the soot growth parameters 
are calculated off-line is the relative case by which different growth models can be incorporated into 
the analysis, and the fact that one is not forced to oversimplify either the chemistry or the soot 
growth mechanism as is presently the case with transport-based approaches. The Lindstedt growth 
parameters used are tabulated below. 

Table 3 

Rate Constants, in the form A T b e _E / RT , for Soot Formation and 
Consumption Model (units are kg, m, s, kmol, kcal, and K) 


Rate Constant 

A 

b 

E 

k* 

1.35 x 10 6 

0.0 

41 x 10 3 

k a 

5.00 x 10 3 

0.0 

24 x 10 3 

kiii 

1.78 x 10 4 

0.5 

39 x 10 3 


where the rate constants are associated with inception, surface growth, oxidation and coalescence 
processes in the way described in Figure 15. The model in its original form solves two coupled 
conservation equations, one for the soot mass fraction, and one for the soot number density. The 
present use of these growth parameters differs from that in Ref. 33 in that they are used here with 
the size class representation, which provides a size distribution. Refs. 23-24 assumed a monodisperse 
size distribution. While the results are generally sensitive to inception rate, they are not sensitive 
to the mass selected for the inception species. A nominal mass of 720 a.m.u. has been assumed 
for the inception species. Although the model jet fuel kinetics scheme includes steps leading to 
benzene, only the calculated acetylene concentration is used for purposes of calculating inception 
rates in the Lindstedt model. 

With just two Lindstedt model dependent variables, transport or balance equations for them 
could have been readily incorporated into the turbulent flow solver. (Examples of this can be 
seen in certain papers contained in Ref. 12). However, this is done on a "monodisperse” basis, 
that is, there would only be a single particle size at each location. Using it in conjunction with 
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Figure 15. The Lindstedt global soot formation model. 

(Fairweather, Jones and Lindstedt, Comb. & Flame, 89, 1992) 



the stretched flamelet and size class growth algorithms has the advantage that information about 
particle size distribution is obtained. Also, it is an open question whether monodisperse models 
can properly treat high pressure growth effects associated with different size particles belonging 
to different Knudsen Number regimes. It is felt that the size class representation used here can 
properly treat such effects. Ultimately, the Lindstedt model will probably be basedon a small mass 
PAH inception scheme; stretched flamelet theory, with its ease of incorporating complex chemistry, 
can readily handle this, while transport-based approaches will find it all but impossible. 

Sooting flamelet calculations have been carried out for the conditions given in Table 4 with 
strain rate as an external parameter. 

Table 4 

• 10.5 ATM 

• T / 4=917K (1190F) 

• Tf=478K (400F) 

• Strain rate 2750 - 50 sec -1 

Calculated growth and oxidizing species in the model jet fuel flamelet simulation have been 
shown in Figure 2 for a specific value of strain rate. The temperature profile and volume fraction 
for this case are shown in Figure 16. The soot is seen to be incepted on the fuel side of the flame, 
and grows as it is transported toward the stagnation plane by convection and thermophoresis. The 
average soot particle size profile, reflecting the effects of surface growth and coalescence, is shown in 
Figure 17. Average particle sizes are much larger than those found at atmospheric pressure because 
of the high pressure. The calculation yields additional informatiotn that is not shown here. The 
particle size distribution is calculated at each point, for example, and could be provided. 

The calculated soot volume fraction profiles in mixture fraction space with strain rate as an 
external parameter are shown in Figure 18. As seen, the soot profiles have the form of a similarity 
solution of the form f v (z, a) = h(a)g(z). Numerical differentiation of the profiles to give the scalar 
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Figure 16. Soot volume fraction and temperature profiles in model jet fuel flamelet. 



Figure 17. Average particle diameter in model jet fuel flamelet. 


Log 10 (Soot Volume Fraction) 


Fuel: 15% Trimethylbenzene, 85% Decane 
Lindstedt Model; P = 10.5 atm 



Figure 18. Strain rate dependence of soot formation in model jet fuel flamelets. 




dissipation also reveals a similarity form as discussed in the preceding section on development of 
the joint mixture fraction/scalar dissipation pdf. The flameiet profiles contain all the information 
needed to formulate the joint pdf as shown in the Subtask F discussion. 

While setting up the counterflow code to calculate them the first time is a cumbersome task, 
it is likely that this can be highly automated in the future. Work is in progress on creation of 
“continuation” versions of the opposed jet code, in which the range of strain rates, pressures, fuel 
composition and fuel/air temperatures is swept out continuously (Ref. 25). For pressure, this 
is likely to yield fairly simple scaling laws that can be used to correct parameters obtained at 
a reference pressure. Changing jet fuel kinetics mechanisms also is not a difficult task; utilities 
exist to convert restart files from one mechanism to another. Most of the calculational investment 
in the stretched flameiet approach occurs in the library setup and generation; the postprocessing 
calculation of soot loadings is relatively simple and efficient. 

The predicted pressure dependence of the volume fraction profile for a representative strain 
rate is shown in Figure 19. For low to moderate pressures, soot volume fraction is about quadratic 
in pressure. There are many factors contributing to this dependence. As pressure increases, more 
complete reaction leads to higher overall temperatures in the flameiet, leading in turn to higher 
inception and surface growth through the Arrhenius factors of Table 3. Increased acetylene concen- 
tration also leads to increased inception and surface growth. The result is similar to that observed 
in pre-mixed and co-flow diffusion flames (Refs. 26-29). At much higher pressures, the expecta- 
tion would be that depletion of the gaseous carbon pool would result in a pressure dependence 
approaching a linear relationship. 

Figure 20 shows a comparison of calculations performed using the Lindstedt model with those 
performed using a provisional model from Ref. 22, which has inception linked to calculated benzene 
and phenyl concentrations. The latter model is not as well validated as the Lindstedt model, but 
the differences shown are probably representative of the uncertainties in soot kinetic models at this 
time. Use of another soot growth or jet fuel kinetics model would merely require that the curve fits 
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Figure 19. Pressure dependence of soot formation in model jet fuel flamelets. 



Mixture fraction 


Figure 20. Comparison of soot models in jet fuel flamelets. 


of the shape and similarity parameters described in the Subtask F discussion be redone, a relatively 
simple task. 

COUPLED FLOW, SOOT GROWTH, and RADIATION PROGRAM 

The sooting flamelet and radiation algorithms have been coupled to the TEACH code on a post- 
processing basis. The flow calculation is completed, and the converged parameters are supplied as 
needed input for the joint pdf and radiation calculations. This assumes that the radiation represents 
a small fraction of the total flow enthalpy release. Should the contrary be true in some circumstances, 
it would be possible in principle to perform an iterative calculation in which the flux divergence 
is supplied as a sink term to the energy equation for the next iteration on the flow code, and so 
on. While the TEACH code has been used for demonstration purposes, any other flow code whose 
output can be arranged to give the mean mixture fraction, its variance, and the scalar dissipation 
could be used. 

A flow diagram showing the calculational procedure is shown in Figure 21. The contours 
of the axisymmetric combustor are first supplied to the program TRFN2D, which calculates the 
radiation grid using transfinite interpolation based on the desired number of grid points in the r 
and z directions. A program named READFLOW reads the output of the TEACH program; the 
important quantities for the soot growth and radiation calculation are the mean mixture fraction, 
its variance, and the scalar dissipation at each node point. The program also reads in the means 
of scalars like temperature, density, and the species concentrations, since this version of TEACH 
calculates these, but the single scalar averaging needed for these quantities can also be done in 
the radiation code. The output of READFLOW is then mapped onto the radiative mesh using 
RADMAP, which uses bilinear interpolation. RADMAP also incorporates information on wall 
temperatures and emissivities into its output file. The radiation grid, the output of RADMAP, and 
a flamelet library for single scalar pdf calculations are then input to the discrete transfer program, 
RADCALC. The curve-fits of the sooting flamelet library calculations (Subtask F) are incorporated 
into a subroutine JOINTPDF, which is linked to RADCALC. The options available in RADCALC 
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Figure 21 . Radiation code flow diagram. 
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are a radiation calculation based on time-averaged properties, and a calculation using the turbulent 
radiation algorithm which has been discussed. When the time averaged property option is selected, 
there is also the option of doing narrowband calculations using RADCAL (Ref. 20). The input also 
includes the number of rays to be used in the discrete transfer analysis, in the manner described 
in Refs. 18-19. The output of the program consists of the net radiative fluxes on all walls, and the 
volumetric radiative dissipation rate. Output is in KW/M 2 and KW/M 3 , respectively. Nominally 
32 rays per point per quadrant are used in the discrete transfer analysis. 

Sample calculations have been carried out for the decane-fuelled, bluff-body dump combustor 
configuration shown previously in Figure 3. The geometry is of a type in use at Wright-Patterson 
Air Force Base (Ref. 30). In the base configuration, fuel is injected at the middle of a center-body 
face through a .96 cm diameter tube; the overall diameter of the centerbody is 14 cm, and the 
airstream is located from radial location 7 cm to the outer wall whose radius is 12.7 cm (5 inches). 
The pressure, equivalence ratio, and fuel/air stream temperatures are representative of conditions 
of practical interest. However, the combustor is of a type intended for diagnostic studies and is 
relatively slow mixing; the calculations to follow are intended only as numerical exercises, and are 
not meant to represent a combustor of commercial interest. Experiments in such a simple geometry 
would be ideal for model validation. 

The average temperature distribution in the model calculation is shown in Figure 22. The 
relative slowness of the mixing is indicated by the fact that the average temperature level is still 
rising at two meters. These average temperatures were calculated using a single scalar, mixture 
fraction pdf, since temperature is only slowly varying with strain rate, and there is no need to use 
the joint pdf. Corresponding mean mixture fraction and scalar dissipation distributions are shown 
in Figure 23. The scalar dissipation, in units sec^” 1 ^, is seen to die away relatively rapidly with 
distance from the fuel injector. These quantities, together with the mixture fraction variance, are 
needed for the single scalar and joint pdf averaging algorithms. Application of the joint pdf results 
in the calculated average soot distribution shown in Figure 24. The soot peaks on the centerline 
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Figure 22. Temperature distribution in dump combustor simulation. 
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near the fuel injector at values somewhat in excess of 10^“ and dies away over the course of the 
two meter length. Peak soot volume fractions in ten atmosphere, kerosene-fuelled combustors are 
known to be of this magnitude (Ref. 31); while much more precise theory-experiment validation is 
necessary before any conclusions can be drawn, this is preliminary encouragment that these first 
calculations of soot growth in high pressure combustors are of the right order. (Soot loadings are 
sometimes given in units of grams per cubic meter. The predicted peak soot volume fractions in 
this simulation are on the order of 40 g/m s . There is experimental evidence that primary zone soot 
concentrations at ten atmospheres are indeed about this level (Ref. 31)). 

Averaging out the strain rate in the joint pdf gives the pdf in terms of mixture fraction, 
illustrated at the point of maximum soot in Figure 25. Whereas in an individual flamelet the 
volume fraction maximum occurs at a mixture fraction of about .65, near the stagnation plane, the 
joint pdf maximum occurs at lower mixture fractions around .4. This reflects the effects of flamelet 
stretching by the turbulent flow. 

Net radiative fluxes on the cylindrical wall are shown in Figures 26-28. The pure soot, pure 
gas, and gas plus soot on a time-averaged basis are shown in Figure 26. Soot radiation is seen to 
dominate the gas radiation for soot levels of this magnitude. The gas radiation is dominated by 
the 4.3 micron band of CO 2 . The good agreement of the present radiation model with spectrally- 
integrated, narrowband calculations is shown in Figure 27, again on a time-averaged basis. The 
present wideband-based model of Equations 25, 31 and 32 is much more efficient than the narrow- 
band calculations, and gives agreement that is entirely satisfactory. Comparison to a prediction 
using the turbulent radiation algorithm shows that in this case, turbulence is predicted to result in 
an enhancement of somewhat less than 20% relative to calculations based on time-averaged prop- 
erties (Figure 28). Because the double integration involved in the joint pdf is at present quite 
time-consuming, the turbulent radiation algorithm as it relates to soot is in the code in an approx- 
imate form. The temperature is allowed to fluctuate according to the single scalar pdf, using the 
average volume fraction calculated with the joint pdf. A priority item for future work must be to 
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Figure 25. Joint mixture fraction/scalar dissipation PDF in model jet fuel/ 
dump combustor simulation. 




Figure 26. Soot and gas radiation fluxes on dump combustor wall. 
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Figure 28. Turbulence enhancement of radiation on dump combustor walls. 



find a more efficient double integration technique. 


The radiative dissipation profile on a time-averaged basis is shown in Figure 29, corresponding 
to the contributions of both gas and soot. The dissipation rate (net radiative emission rate) tends to 
follow the temperature distribution to an exaggerated extent because of the sensitivity of the Planck 
function to temperature (Compare Figure 22). Because of the time-consuming nature of the joint 
pdf calculation, radiative dissipation is presently calculated only on a time-averaged property basis; 
only the fluxes on the walls are calculated using the turbulent radiation algorithm, and not the 
internal fields needed for the dissipation calculation. The radiative source term could be supplied 
as an energy sink term to the flow code energy equation for an iterative calculation to see whether 
radiation significantly depresses average temperatures. The dissipation rate is the local emission 
rate minus the rate of absorption of radiation from other parts of the combustor; note in Figure 
29 a region of cold soot on the centerline near the fuel injector in which there is net absorption of 
radiation. 

Extensive parametric variations have not yet been carried out with this model. Changing the 
fuel and air flow rates by the same factors is found to result in a mixture fraction /mixture fraction 
variance profile that is sensibly unchanged in spatial coordinates. The peak scalar dissipation does 
scale with velocity in the manner expected, but its influence is confined to a relatively small region 
near the fuel nozzle. The result is a soot profile that is relatively insensitive to the velocities of 
the streams (keeping equivalence ratio constant). The radiative wall fluxes are similarly insensitive, 
implying that the fraction of the total enthalpy converted to radiation is inversely proportional 
to velocity at constant equivalence ratio. The predicted pressure dependence of the soot volume 
fraction would approximate that of the flameiets; that is, it would be about quadratic for pressures 
in the vicinity of ten atmospheres or below, and would be expected to make a transition from 
quadratic to linear at much higher pressures. 

A sample, provisional, application of this soot formation theory to the RBQQ sector rig is 
shown in Figure 30. A three-dimensional CFD flowfield simulation was obtained courtesy of CFD 
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Research, Inc. From the provided profiles of mean mixture fraction, its variance, and the scalar 
dissipation, values were extracted along the center of the sector rig, and soot loadings calculated 
using the joint PDF algorithm. As seen, peak soot volume fractions in the rich zone approach 10 5 . 
The soot is seen to be very effectively oxidized in the quench zone. 

CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK 

These investigations have lead to a number of noteworthy technical achievements. Among these 
is the first calculation of soot formation in jet fuel based on a complex chemical mechanism. The 
associated soot growth calculations are the first to incorporate continuum effects in the particle 
kinetics scheme in order to treat high pressure growth. Further, this work marks the first inclusion 
of efficient turbulence-radiation interactions algorithms into a radiation code. The inclusion of non- 
linear effects in the k-e flow code and the formulation of a joint mixture fraction/scalar dissipation 
pdf for turbulent soot formation from the stretched flamelet calculations are also significant technical 
aspects. In the model jet fuel, predicted soot levels appear to be of the right order when compared 
to primary zone data at elevated pressure. However, these data are limited; before this analysis 
can be applied with confidence as a design tool, future work directed towards model validation and 

enhancement should be undertaken. 

Model validation studies would involve comparisons with turbulent jet data, starting with 
simpler fuels like ethylene and propane. In terms of other basic experimental data, there is a strong 
need for high pressure, sooting opposed jet experiments, starting with ethylene, and progressing to 
more complex fuels. Well-diagnosed model combustor experiments providing soot volume fractions 
and radiative fluxes also would be invaluable. 

Certain model enhancements could be undertaken. Among these would be sensitivity analyses 
and simplification of the jet fuel kinetics mechanism. Simplified correlations and scaling relation- 
ships linking soot levels to pressure, strain rate, and fuel/air temperatures can be developed to 
obviate the need to generate new flamelet libraries. As new information on soot inception and sur- 
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face growth becomes available, these would be evaluate for impact on jet fuel soot formation. The 
basic stretched flamelet theory could be enhanced by consideration of partial premixing, flamelet- 
flamelet interactions, and non-adiabatic loss. 

While the TEACH code was employed in these demonstration calculations, any flow code 
that can be configured to provide the required joint pdf parameters could be employed. Thus, 
for example, extensions to unstructured grid codes such as CORSAIR would not be a complicated 
matter. The discrete transfer radiation algorithm used here has considerable geometric flexibility. 
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NOMENCLATURE 


a beta density parameter; strain rate in opposed jet flame 

A band absorptance 

A' band absorptance derivative 

b beta density parameter 

ci, C2 Planck function constants 

c a factor in soot absorption coefficient 

D soot diameter 

Ei exponential integral 

f v soot volume fraction 

H height above burner surface 

I radiative intensity 

lb Planck function 

k absorption coefficient 

p probability density for mixture fraction 

Ro fuel tube radius 

s optical pathlength 

T gas temperature 

z fuel mixture fraction 

z ' fluctuation in mixture fraction 
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Greek symbols 

a integrated band intensity 

X scalar dissipation 

A Length of locally homogeneous portion of inhomogeneous path 

e emissivity 

7 parameter in beta density 

7E Euler-Mascheroni constant 

T Gamma function 

A thermal radiation wavelength 

rji probability density parameter 

A uj bandwidth 

£ band intensity path integral 

p radiating gas density 

(jj frequency 

u>(°) band center frequency 

Subscripts 

i i_th molecular resonance 

w wall value 
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APPENDIX A 


Reaction Mechanism Rate Coefficients in the Form kf = AT 0 exp(—Eo/RT). 
Units are moles, cubic centimeters, seconds, Kelvins, and calories/mole. 



REACTION 


fi 

_E_ 

1. 

H+02=0+0H 

5.10E16 

-0.820 

16510. 

2. 

H2+0=H+0H 

1.80E10 

1.0 

8830. 

3. 

H2+0H=H20+H 

1.20E09 

1.3 

3630. 

4. 

0H+0H=H20+0 

6.0E08 

1.3 

0 . 

5. 

H+0H+M=H20+M (M=AR) 
H2O/20./ 

7.50E23 

-2.6 

0 . 

6. 

02+M=0+0+M 

1.90E11 

0.5 

95560. 

7. 

H+H+M=H2+M (M=AR) 
H20/0.0/H2/0.0/C02/0.0/ 

1.0E18 

-1.0 

0 . 

8. 

H+H+H2=H2+H2 

9.20E16 

-0.6 

0 . 

9. 

H+H+H20=H2+H20 

6.00E19 

-1.250 

0 . 

10. 

H+H+C02=H2+C02 

5.49E20 

-2.0 

0 . 

11. 

H2+02=0H+0H 

1.70E13 

0.0 

47780. 

12. 

H+02+M=H02+M (M=AR) 
H2O/21./CO2/5./H2/3.3/CO/2./O2/0./N2/0./ 

2.10E18 

-1.0 

0 . 

13. 

H+02+02=H02+02 

6.70E19 

-1.420 

0 . 

14. 

H+02+N2=H02+N2 

6.70E19 

-1.420 

0 . 

15. 

H02+H=H2+02 

2.50E13 

0 . 

700. 

16. 

H02+H=0H+0H 

2.50E14 

0 . 

1900. 

17. 

H02+0=0H+02 

4.80E13 

0 . 

1000. 

18. 

H02+0H=H20+02 

5.00E13 

0 . 

1000. 

19. 

H02+H02=H202+02 

2.00E12 

0 . 

0 . 

20. 

H202+M=OH+OH+M 

1.20E17 

0 . 

45500. 

21. 

H202+H=H02+H2 

1.70E12 

0.0 

3750. 

22. 

H202+0H=H20+H02 

1.00E13 

0.0 

1800. 

23. 

C0+0-(-M=C02+M 

3.20E13 

0.0 

-4200. 

24. 

C0+02=C02+0 

2.50E12 

0 . 

47700. 

25. 

C0+0H=C02+H 

1.50E07 

1.3 

-760. 

26. 

C0+H02=C02+0H 

5.80E13 

0 . 

22930. 

27. 

CH4+M=CH3+H+M (M=AR) 
H20/5./ 

1.00E17 

0 . 

88000. 

28. 

CH4+H=CH3+H2 

2.20E04 

3.0 

8750. 

29. 

CH4+0=CH3+0H 

1.20E07 

2.080 

7630. 

30. 

CH4+0H=CH3+H20 

3.50E03 

3.080 

2000. 

31. 

CH4+CH2=CH3+CH3 

1.30E13 

0 . 

9500. 

32. 

CH3+M=CH2+H+M 

1.90E16 

0 . 

91600. 

33. 

CH3+CH3=C2H6 

1.60E13 

0.0 

-306. 

34. 

CH3+CH3=C2H4+H2 

2.10E14 

0.0 

19200. 

35. 

CH3+CH2=C2H4+H 

3.00E13 

0.0 

0 . 
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36. 

CH3+H=CH2+H2 

9.00E13 

0 . 

15100. 

37. 

CH3+0=CH20+H 

6.80E13 

0 . 

0 . 

38. 

CH3+0=CH2+0H 

5.00E13 

0 . 

12000. 

39. 

CH3+0H=CH2+H20 

1.50E13 

0 . 

5000. 

40. 

CH3+OH=CH20+H2 

1.00E12 

0 . 

0 . 

41. 

CH3+02=CH20+0H 

5.20E13 

0 . 

34570. 

42. 

CH3+02=CH30+0 

7.00E12 

0 . 

25650. 

43. 

CH30+M=CH20+H+M 

1.00E14 

0 . 

25000. 

44. 

CH30+H=CH20+H2 

2.00E13 

0.0 

0 . 

45. 

CH30+0=CH20+0H 

1.00E13 

0 . 

0 . 

46. 

CH30+0H=CH20+H20 

1.00E13 

0 . 

0 . 

47. 

CH30+02=CH20+H02 

6.30E10 

0 . 

2600. 

48. 

CH20+M=HC0+H-fM 

3.31E16 

0.0 

81000. 

49. 

CH20+H=HC0+H2 

2.20E08 

1.770 

10500. 

50. 

CH20+0=HC0+0H 

1.80E13 

0.0 

3080. 

51. 

CH20+0H=HC0+H20 

3.40E09 

1.180 

-447. 

52. 

HCO+M=CO+H+M 

1.60E14 

0 . 

14700. 

53. 

HCO+H=CO+H2 

4.00E13 

0.0 

0.0 

54. 

HC0+0=C0+0H 

3.00E13 

0.0 

0.0 

55. 

HC0+0=C02+H 

3.00E13 

0.0 

0.0 

56. 

HC0+0H=C0+H20 

5.00E12 

0.0 

0.0 

57. 

HC0+02=C0+H02 

3.30E13 

-0.4 

0.0 

58. 

CH2+H=CH+H2 

7.30E17 

-1.560 

0 . 

59. 

CH2+0=C0+H+H 

3.00E13 

0.0 

0 . 

60. 

CH2+0=C0+H2 

5.00E13 

0.0 

0 . 

61. 

CH2+0=CH+0H 

5.00E13 

0.0 

12000. 

62. 

CH2+0H=CH20+H 

3.00E13 

0.0 

0 . 

63. 

CH2+OH=CH+H20 

4.50E13 

0.0 

3000. 

64. 

CH2+02=C02+H+H 

1.60E12 

0.0 

1000. 

65. 

CH2+02=C02+H2 

6.90E11 

0.0 

500. 

66. 

CH2+02=C0+H20 

1.90E10 

0.0 

-1000. 

67. 

CH2+02=CO+OH+H 

8.60E10 

0.0 

-500. 

68. 

CH2+02=HCO+OH 

4.30E10 

0.0 

-500. 

69. 

CH2+02=CH20+0 

2.00E13 

0.0 

9000. 

70. 

CH2+C02=C0+CH20 

1.10E11 

0.0 

1000. 

71. 

CH+0=C0+H 

5.70E13 

0.0 

0 . 

72. 

CH+OH=HCO+H 

3.00E13 

0.0 

0 . 

73. 

CH+02=HC0+0 

3.30E13 

0.0 

0 . 

74. 

CH+C02=HC0+C0 

3.40E12 

0.0 

690. 

75. 

C2H6+H-C2H4+H+H2 

5.40E02 

3.5 

5200. 

76. 

C2H6+0H-C2H4+H+H20 

8.70E09 

1.050 

1810. 

77. 

C2H6+CH3-C2H4+H+CH4 

5.50E-01 

4.0 

8280. 

78. 

C2H4+M=C2H2+H2+M 

2.60E17 

0.0 

79350. 

79. 

C2H4+M=C2H3+H+M 

2.60E17 

0.0 

96600. 

80. 

C2H4+H=C2H3+H2 

1.10E14 

0.0 

8500. 

81. 

C2H4+0H=C2H3+H20 

4.80E12 

0.0 

1230. 

82. 

C2H4+0H=CH20+CH3 

2.00E12 

0.0 

960. 
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83. 

C2H3+M=C2H2+H+M 

8.00E14 

0.0 

31500. 

84. 

C2H3+H=C2H2+H2 

4.00E13 

0.0 

0 . 

85. 

C2H3+02=HC0+CH20 

4.00E12 

0.0 

-250. 

86. 

C2H2+M=C2H+H+M 

4.20E16 

0.0 

107000. 

87. 

C2H+H2=C2H2+H 

4.10E05 

2.390 

860. 

88. 

C2H2+0=CH2+C0 

2.20E10 

1.0 

2580. 

89. 

C2H2+OH=CH2CO+H 

3.20E11 

0.0 

200. 

90. 

C2H2+0H=C2H+H20 

6.00E12 

0.0 

7000. 

91. 

CH2CO+M=CH2+CO+M 

3.60E15 

0.0 

59300. 

92. 

CH2CO+H=CH3+CO 

1.10E13 

0.0 

3430. 

93. 

CH2C0+0=CH20+C0 

2.0E13 

0.0 

0 . 

94. 

CH2C0+0H=CH20+HC0 

2.80E13 

0.0 

0.0 

95. 

C2H+0=CH+CO 

5.00E13 

0.0 

0 . 

96. 

C2H+02=CO+HCO 

2.40E12 

0.0 

0 . 

97. 

C10H22+H=C10H21+H2 

(WARNATZ) 

3.0E14 

0.0 

8457. 

98. 

C10H22+O=C10H21+OH 

(WARNATZ) 

1.0E14 

0.0 

4539. 

99. 

C10H22+OH=C10H21+H2O 

(WARNATZ) 

1.0E13 

0.0 

883. 

100. 

C10H22+CH3=C10H21+CH4 

(PITZ) 

2.0E11 

0.0 

9500. 

101. 

C10H22+C2H3=C10H21+C2H4) 

(PITZ) 

3.0E11 

0.0 

1800. 

102. 

C10H22+02=C10H21+H02 

2.51E13 

0.0 

49000. 

103. 

C10H21-C5H11+C5H10 

(WARNATZ) 

2.5E13 

0.0 

28715 

104. 

C5H11-2C2H4+CH3 

2.5E13 

0.0 

28715 

105. 

C5H10+H=C2H4+C3H5+H2 

(WESTBROOK) 

5.0E13 

0.0 

3900. 

106. 

C5H 10+O=C2H4+C3H5+OH 
(WESTBROOK) 

5.0E13 

0.0 

3900. 

107. 

C5H 10+OH=C2H4+C3H5+H2O 
(WESTBROOK) 

5.0E13 

0.0 

3900. 

108. 

C5H10+CH3=C2H4+C3H5+CH4 

(WESTBROOK) 

5.0E13 

0.0 

3900. 

109. 

C5H10+C2H3=2C2H4+C3H5 

(WESTBROOK) 

5.0E13 

0.0 

3900. 

110. 

C3H5=C3H4+H 

3.98E13 

0.0 

70000. 

111. 

C3H5+02=C3H4+H02 

6.03E11 

0.0 

10000. 

112. 

C3H4+0=CH20+C2H2 

1.0E12 

0.0 

0.0 

113. 

C3H4+0H=CH20+C2H3 

1.0E12 

0.0 

0.0 

114. 

C3H4+H=CH3+C2H2 

2.0E13 

0.0 

2411 

115. 

C3H4+0=C0+C2H4 

1.4E13 

0.0 

2103 
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116. 

C7H8+02-C6H6+CH20+0 

(GUERRET) 

1.0E18 

0.0 

50000. 

117. 

C8H10+02-C7H8+CH20+0 

(GUERRET) 

1.0E18 

0.0 

50000. 

118. 

C9H12+02-C8H10+CH20+0 

(GUERRET) 

1.0E16 

0.0 

50000. 

119. 

C6H6=PHENYL+H 
(JACKSON AND LAURENDAU) 

5.0E15 

0.0 

108000. 

120. 

C6H6+H=PHENYL+H2 
(JACKSON AND LAURENDAU) 

3.0E12 

0.0 

8100. 

121. 

C6H6+0=PHENYL+0H 

2.8E13 

0.0 

4910. 

122. 

C6H6+0H=PHENYL+H20 

2.1E13 

0.0 

4570. 

123. 

PHENYL+02-2CO+C2H2+C2H3 

1.1E15 

0.0 

24000. 
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APPENDIX B 


The variance of Xst, VAR, is given by the expression 

VAR = (x. 2 t > - <x.t> 2 

where (x 2 t ) and (x»t) are the first and second moments of x.t, respectively. 


<X.t) 2 = f xft P(x.t)dx.t 

Jo 

1*00 

(x»t> = / x*t P(x*t)dx.t 

JO 


Substituting the expression for the marginal density, 


and 


(x.t ) 2 = f xlt f P(x.t,z)dz dx,t 

Jo JO 

(X»t) = f X>t f P(x«t,z)dz dx.t 
Jo Jo 


(B- 1 ) 


(B- 2) 


(B~ 3) 


Substituting x/f for x»t and changing the integration variable yields the following expressions 
for for the first two moments, 


{X.0 2 = (X 2 )^ 1 ^P(z)dz (fl-4) 

(x*t) = (x) / fTT P(z)dz 
Jo H Z J 

(x 2 ) and (x) "e found through the analytical expression for the moments of the lognormal distri- 
bution, (Ref. ) 

(X*) = exp(2ji + ff 2 ) • (exp <r 2 - 1) [B - 5) 

(x) = exp(/z + i<r 2 ) 

Based on experimental data, a is assumed equal to unity, (x) is computed in the mean flow 
calculation; thus, n is known. The integrals in z space can be evaluated analytically to compute 
VAR at each point. However, for purposes of estimating the interval of integration for Xit, the 
following approximations were made, 

Jo f Hz) P(z)dz = ( f*(^ “ f*<i) (B ~ 6) 
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Hence, 


VAR = 


and the standard deviation of Xst,SD, is 


SD = VAR 1/2 


(B-7) 


It was found an integration interval from x»t = 0 to Xst +4SD was sufficient. 
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